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ABSTRACT 

We measure long (2200-4000 A) and short (1450-2200 A) wavelength spectral slopes a [F^ cx v"") 
for quasar spectra from the Sloan Digital Sky Survey. The long and short wavelength slopes are 
computed from 3646 and 2706 quasars with redshifts in the z=0.76-1.26 and z=1.67-2.07 ranges, 
respectively. We calculate mean slopes after binning the data by monochromatic luminosity at 2200 
A and virial mass estimates based on measurements of the Mg II line width and 3000 A continuum 
luminosity. We find little evidence for mass dependent variations in the mean slopes, but a significant 
luminosity dependent trend in the near UV spectral slopes is observed with larger (bluer) slopes at 
higher luminosities. The far UV slopes show no clear variation with luminosity and are generally lower 
(redder) than the near UV slopes at comparable luminosities, suggesting a slightly concave quasar 
continuum shape. We compare these results with Monte Carlo distributions of slopes computed from 
models of thin accretion disks, accounting for uncertainties in the mass estimates. The model slopes 
produce mass dependent trends which are larger than observed, though this conclusion is sensitive 
to the assumed uncertainties in the mass estimates. The model slopes are also generally bluer than 
observed, and we argue that reddening by dust intrinsic to the source or host galaxy may account for 
much of the discrepancy. 

Subject headings: accretion, accretion disks — black hole physics — galaxies: active — galaxies: 
quasars: general 



1. INTRODUCTION 

The "bare" thin accretion disk has been invoked ubiq- 
uitously to explain both general and detailed properties 
of a wide variety of accreting systems. In systems which 
are believed to harbor black holes the results have been 
mixed. The timing and spectral variations of Galactic 
black hole candidates show a range of behavior, much of 
which cannot be accommodated by a simple thin accre- 
tion disk. However, some of these sources enter a thermal 
state, where the spectral energy density (SED) is dom- 
inated by emission which can be well fit with a simple 
multitemperature blackbody model (Mitsuda et al. 1984) 
or more elaborate variations (e.g. Li et al. 2005; Davis & 
Hubeny 2006). 

The thin disk model may also explain a portion of 
the emission in Active Galactic Nuclei (AGNs). These 
sources are almost certainly powered by an accretion 
flows onto a central "super-massive" black hole (Krolik 
1999). As a result, it is often supposed that the broad 
UV peak ( "big blue bump" ) of the SED is emission from 
a thin accretion disk. However, this interpretation faces 
a number of challenges when detailed comparisons are 
made between models and the data (see e.g. Koratkar 
& Blaes 1999). It may be possible to resolve some of 
these discrepancies by modifying the model and consid- 
ering additional processes (e.g. dust reddening, irradia- 
tion, inhomogeneities, and multi-phase flows), but it is 
then pertinent to ask what, if any, predictive power is 
provided by the thin disk model itself. 

Models of thin accretion disks predict that the effective 
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temperature at the inner-most radius Tin of the disk will 
be proportional to the one-fourth power of accretion rate 
M and inversely proportional to one-half power of black 
hole mass M [T^^ cx M/A'P). This scaling follows rather 
generally from the assumptions that gravitational bind- 
ing energy is radiated locally and that the relevant length 
scale of the emitting region is the gravitational radius 
Rg = GM/c^ of the black hole. As we shall discuss, this 
relation accounts for much (but not all) of the spectral 
dependence on M and M since Tin roughly determines 
the photon energy where the SED peaks (~ A:BTi„). 

In fact, one of the successes of this scenario is that 
it approximately predicts the position of the continuum 
peak for both super-massive black holes in AGN and the 
~ IOA/0 black holes in Galactic X-ray binaries, assum- 
ing that both sources accrete mass at similar fraction 
of the Eddington limit. Also, the approximate relation 
(X L (X M (where L is the bolometric luminosity) can 
be inferred from the spectral evolution of several black 
hole binaries in thermal state (see e.g. Gierlinski & Done 
2004). Unfortunately the narrow range of dynamically 
inferred masses, uncertainties in the estimates, and small 
sample of sources complicate efforts to simultaneously 
and independently constrain the M dependence of the 
SED in these systems. 

In many respects, AGNs offer greater promise for test- 
ing the M and L (or M) dependence of thin disk spec- 
tral models since there is a much larger sample available 
which spans a wider range of Af and L. However, there 
are also a number of additional challenges. First, mass 
estimates in AGN are generally more uncertain than in 
black hole candidates, and the most reliable methods can 
only be applied in a small fraction of sources. Also, the 
SED peaks in the far UV so it is challenging to get broad- 
band coverage at low redshifts, except for a relatively 
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small number of bright sources. Finally, the large num- 
ber of emission lines which characterize most AGN may 
prevent a robust determination of underlying continuum 
emission. 

Despite these challenges, comparing the predicted and 
observed spectral evolution as a function of M and L is 
the principle aim of this work. To accomplish this, we 
measure UV spectral slopes (a where (x u"") for sev- 
eral thousand quasars from the Sloan Digital Sky Survey 
(SDSS). A wide range of a is inferred, allowing us to ex- 
amine the extent to which a correlates with the M and 
L, the parameters which determine the model SEDs (see 
e.g. Sun & Malkan 1989; Laor & Netzer 1989; Hubeny 
et al. 2000). We compare the observed slopes with cal- 
culations based on the relativistic, fully non-LTE models 
of Hubeny et al. (2000). In addition to allowing us to 
test the predictions of the thin disk model, a correlation 
(or lack thereof) may provide clues and place important 
constraints on other processes which may be important 
for determining the continuum SED. 

The plan of this work is as follows. In order to test 
the predictions of the thin disk model we construct a 
large table of artificial SEDs for direct comparison with 
data. The methods used to construct these models are 
summarized in §2. In §3.1 we briefly review the mass 
estimation methods employed in this work as our con- 
clusions are sensitive to the reliability of these estimates. 
In §3.2 we compare the model SEDs with a small sample 
of well-observed, relatively nearby AGNs with simulta- 
neous optical to UV spectra (Shang et al. 2005). In §4 we 
present slope measurements for a large sample of SDSS 
QSOs. In §5 we calculate slopes from our spectral models 
and generate Monte Carlo realizations of the slope distri- 
butions for comparison with the data. In §6 we discuss 
the possible origins of discrepancies between the models 
and observations, considering additional processes, not 
accounted for by the models, which might be important. 
In §6.4, we compare our results with previous work, par- 
ticularly that of Bonning et al. (2006) which is the most 
similar to our current efforts. We summarize our conclu- 
sions in §6.5. Throughout this work, we use the follow- 
ing cosmology: Hq = 70 kms~^ Mpc~^, = 0.3, and 
l^A = 0.7. 

2. SPECTRAL MODELS 

Our artificial SEDs are based on time-independent 

models of thin, a-disks (Shakura & Sunyaev 1973). Here, 
a refers to an assumed constant of proportionality be- 
tween the accretion stress and total pressure. Therefore, 
it is a dimensionless parameter and generally assumed to 
be less than unity. (Hereafter, we will use ass to diff'er- 
entiate this quantity and the spectral slope a, but we will 
continue to refer to the model simply as an a-disk.) We 
generate artificial SEDs using AGNSPEC, an interpola- 
tion scheme which is identical to that described in Davis 
& Hubeny (2006). Except for the interpolation scheme, 
the models used here are equivalent to those presented 
in Hubeny et al. (2000). We only summarize the most 
relevant features and the reader is referred to this work 
(and references therein) for a more detailed discussion. 

The artificial SEDs are based on a relativistic, thin, a- 
disk model similar to that of Novikov & Thorne (1973), 
but include some minor corrections (Riffert & Herold 
1995). These SEDs account for the effects of light bend- 



ing and time dilation by calculating the null geodesies 
of the black hole spacetime (KERRTRANS, Agol 1997). 
In addition to M and M, such fully-relativistic mod- 
els require a choice of the black hole spin parameter 
a* = cJ/GM'^, where J is the angular momentum of 
the hole. 

Using these "one-zone" disk models as a basis, we com- 
pute two types of artificial SEDs. We first construct rela- 
tively simple spectra which assume the disk surface emits 
like a blackbody at the local effective temperature (but 
still account for relativistic effects) . A second set of more 
sophisticated models employ TLUSTY (Hubeny & Lanz 
1995) to solve the coupled equations of radiative trans- 
fer and equilibrium, non-LTE vertical structure in the 
disk. In this case the disk surface density is needed to 
determine the structure, requiring a model for angular 
momentum transport. We adopt an a-disk prescription 
with ass = 0.01. At the wavelengths of interest the 
spectra are relatively insensitive to ass for the ranges of 
parameters (a», M, and L) explored here. These mod- 
els include bound-free and free-free opacities of H and 
He and the effects of electron scattering are calculated 
in the Thomson limit. Due to the difficulties involved, 
we do not include opacity from bound-bound transitions, 
although it may have a significant effect on the resulting 
spectra (see e.g. Hubeny & Hubeny 1998). 

The inclusion of non-LTE effects and realistic opacities 
can cause significant shifts from a simple blackbody spec- 
trum. The effects are usually largest at frequencies near 
the bound-free transitions. Since we concentrate on the 
SED at wavelengths longer than 1000 A, the Balmer edge 
will be the most important. Unfortunately, the annuli 
which produce most of the fiux longward of the Balmer 
edge are the most difficult to calculate due the presence 
of ionization zones associate with the transition. In these 
annuli, equilibrium solutions are either unattainable or 
lead to unstable atmospheres with density inversions (see 
§3.3 of Hubeny et al. 2000). Due to these difficulties, we 
simply assume blackbody spectra for annuli with effec- 
tive temperatures below 9,000 K. Based on Fig. 11 of 
Hubeny et al. (2000), this probably leads to a slight un- 
derestimate of flux at 4000 A and possible implications 
are discussed in §6.3. 

3. METHODS 
3.1. Mass Estimates 

Stellar dynamical estimates of black hole masses are 
only available in the local universe. At higher redshifts 
virial mass estimates based on reverberation mapping of 
broad line region (BLR) clouds (see e.g. Peterson et al. 
2004) are generally accepted to be the most reliable (see, 
however, Krolik 2001). This technique uses the fuU- 
width-at-half-maximum (FWHM) or second moment of 
one or more prominent broad emission lines to estimate 
the velocity field wblr of broad line clouds. Reverbera- 
tion mapping then provides a characteristic radius i?BLR 
from which the virial mass Afvir ^ v^^y{Rbw./G can be 
estimated. The precise normalization depends on the 
kinematics of the BLR and cannot, in general, be deter- 
mined reliably for individual sources. A single normal- 
ization for all sources can be obtained by requiring virial 
estimates to lie on the M—a relation (Onken et al. 2004). 

A significant drawback of this method is that it re- 
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TABLE 1 
Slope Comparison 



Source 


Mass'' 


a(1450 - 22OOA) 


a(2200 - 4OOOA) a(1450 - 4000A) 




(1O»M0) 


Data a = = 0.9 


Data = a = 0.9 Data a = a = 0.9 



3C 273 8.86 

PG 0953+414 2.76 

PG 0052+251 3.69 

Ton 951 0.92 

Mrk 509 1.43 



-0.64 


-0.34 


-0.24 


-0.16 


-0.17 


-0.12 


-0.60 


-0.21 


-0.038 


-0.98 


0.00 


0.073 


-1.0 


-0.074 


0.092 



0.08 


-0.16 


-0.11 


-0.27 


-0.12 


-0.080 


-0.28 


-0.28 


-0.14 


-0.57 


-0.14 


-0.074 


-0.19 


-0.24 


-0.096 



-0.21 


-0.23 


-0.16 


-0.23 


-0.14 


-0.095 


-0.41 


-0.25 


-0.097 


-0.74 


-0.08 


-0.014 


-0.52 


-0.17 


-0.019 



^Reverberation mapping estimates from Peterson et al. (2004). 



quires sources to be frequently monitored and can only 
be robustly applied in cases where the lines provide a 
clear response to continuum variations. As a result, 
it has only been used successfully for relatively nearby 
sources. (Wandel et al. 1999; Kaspi et al. 2000, 2005). 
At higher redshifts, empirically calibrated luminosity- 
radius relations are more commonly used (e.g. Vester- 
gaard 2002; Woo & Urry 2002; Vestergaard & Peterson 
2006). For these methods it is generally assumed that 
-Rblr, fx (XLxY where monochromatic luminosity XLx 
is calculated at some particular rest wavelength. This 
is usually chosen to be 5100 A when wblr is estimated 
from the width of H/3. The constant of proportionality 
and exponent S are then fit to minimize the scatter be- 
tween this relation and estimates from a reverberation 
mapped sample. Typically S ~ 0.5 (see e.g. Bentz et al. 
2006) is obtained, consistent with expectations from pho- 
toionization models of the BLR. 

At higher redshifts, H/3 shifts out of the observed band 
and other broad emission lines are used, typically Mg II 
or C IV. The estimates used in §4 were calculated from 
line width measurements obtained by McLure & Dunlop 
(2004) and the reader is referred to their paper for a de- 
tailed discussion. Since we are only considering quasars 
with ;z > 0.7 all estimates utiUze the FWHM of Mg II 
WMg and L3000 the monochromatic luminosity measured 
at 3000 A. Therefore, the mass estimates arc calculated 
using equation (A6) of McLure & Dunlop (2004) 




Using L3000, McLure & Jarvis (2002) obtain 6 = 0.47 
by fitting the -Rblr, ^ ^3000 relation to best match re- 
verberation mapped masses of 34 AGNs (Kaspi et al. 
2000). However, this relation is derived including a num- 
ber of relatively low luminosity AGNs while the z > 0.7 
SDSS quasar sample only includes sources with L3000 > 
10'''* erg s-i. Therefore, McLure & Dmilop (2004) refit 
the i?BLR ^ -^3000 relation only using those sources in the 
reverberation mapped sample with L3000 > 10^^ erg s^^, 
finding 6 = 0.62. Since our sample is a subset of the 
McLure & Dunlop (2004) sample which only includes 
this luminosity range, we also adopt the S = 0.62 value. 

3.2. Comparing Models to Data with Spectral Slopes 

One of the most striking characteristics of most AGNs 
(particularly quasars and type I Seyferts) are the numer- 
ous strong broad, emission lines. A substantial fraction 



of the UV and optical fiux can be emitted in these lines, 
making a precise and robust identification of the contin- 
uum difficult, if not impossible. In some cases the broad 
line region emission appears to be unpolarized and the 
underlying continuum can be identified with polarime- 
try (Kishimoto et al. 2003, 2004), revealing a Balmer 
edge beneath the "small blue bump" . Although promis- 
ing, this technique requires polarimetry observations and 
unpolarized BLR emission, making it unsuitable for our 
purposes. Alternatively, we could attempt to model the 
BLR emission. Photoionization models can qualitatively 
reproduce many aspects of the observed spectra, but de- 
tailed fits to individual SEDs remain a challenge. There- 
fore, we feel the substantial work involved in generating 
such models and fitting the data is best left for future 
efforts. 

Lacking a suitable method for removing the BLR emis- 
sion, we simply adopt "continuum windows" which ap- 
pear to be largely devoid of BLR contamination and fo- 
cus on the flux at these wavelength for comparison with 
the models. Unfortunately, no part of the spectrum is 
completely devoid of line emission, and the best one can 
hope to do is choose regions which are relatively free of 
BLR contamination. Our choices are guided in part by 
the SDSS quasar composite of Vanden Berk et al. (2001). 
By constructing a composite spectrum from thousands 
of quasars, one can construct a high signal-to-noise SED 
which is sensitive to the presence of rather weak emission 
features. We focus on three windows near 1450 A, 2200 
A, and 4000 A. 

In making these choices, we have attempted to bal- 
ance the requirement of little contamination with a desire 
to have as large a sample of quasar spectra as possible 
which cover two of the windows. Of these choices, 2200 
A is probably the most problematic due to the presence 
of Fe II contamination (see e.g. Fig. 6 of Vanden Berk 
et al. 2001). However, there are no preferable options 
between 1450 A and 4000 A. The SDSS spectra cover 
an observed wavelength range from 3800 A to 9200 A, 
roughly a factor 2.4. Therefore, 1450 A and 4000 A can- 
not be simultaneously covered with SDSS, making an 
intermediate window necessary. 

Although we do not have full UV spectral coverage for 
the SDSS quasars, we can gain insight by first compar- 
ing our models with broadband spectra of lower redshift 
AGN. The sample of Shang et al. (2005) is ideally suited 
for this purpose. It combines nearly simultaneous spectra 
for 17 relatively bright, nearby AGN taken with Far Ul- 
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Fig. 1. — Optical to UV SEDs of five AGN from the sample of Shang ct al. (2005). The dashed curves are model SEDs with M 
corresponding to the the reverberation mapping estimates in Table 1 and a* = 0. We choose the inclinations i so that cos i = 0.8 for 
all models except 3C 273. In this case we use cosi = 1 since the source exhibits superluminal motion. The accretion rate is fixed by 
approximately matching the model and data at 4000 A. In Ton 951 and Mrk 509 the bare models are clearly a poor approximation so we 
also consider models (dotted covers) which are reddened using an SMC- like extinction curve. See the text in §3.2 for further discussion. 



traviolet Spectroscopic Explorer, Hubble Space Telescope, 
and the 2.1 m Kitt Peak National Observatory. Shang 
et al. (2005) provide a detailcid analysis of the sample, in- 
cluding comparisons with the spectral models described 
in §2. Here we only provide a brief comparison of five 
sources with reverberation mapping mass estimates (Pe- 
terson et al. 2004). 

As shown in Figure 1, the SEDs of these sources were 
shifted to rest wavelength and corrected for Galactic red- 
dening using the values from Table 1 of Shang et al. 
(2005) which were obtained from NED and based on the 
Schlcgel ct al. (1998) map. The centers of the contin- 
uum windows arc marked with dashed vertical lines. For 
comparison we plot thin disk models for a* = black 
holes as dashed curves. The models can be very roughly 
approximated by power laws except for the breaks near 
the Lyman and Balmer edges at 912 A and 3650 A. The 
broad feature between 2200 A and 4000 A ("small blue 
bump" ) is predominantly a superposition of Balmer con- 



tinuum and a plethora of Fc 11 emission lines produced 
by the BLR, and. therefore, not accounted for bv the 
model SEDs. 

The black hole masses used for the models are equiv- 
alent to the reverberation mapping estimates which are 
listed in the second column of Table 1 (Peterson et al. 
2004). The model M is chosen to approximately match 
the observed flux at 4000 A, assuming an inclination i 
such that cosi = 0.8 for all cases except 3C 273. In the 
case of 3C 273, we used i = since superluminal motion 
is observed in this source. The model bolometric lumi- 
nosities in Eddington units are 0.32, 0.63, 0.063, 0.16, 
and 0.035 for 3C 273, PG 0953-f 414, PG 0052-1-251, Ton 
951, and Mrk 509, respectively. 

Qualitatively, the comparisons provide mixed results. 
They are somewhat better for the higher luminosity 
AGN, but the the model spectra arc clearly too steep in 
the UV for the lowest luminosity sources (Ton 951 and 
Mrk 509). PG 0953-h414 and PG 0052-1-251 yield the 
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best results in 1450-4000 A band, although they respec- 
tively overprcdict and undcrpredict the emission below 
1000 A. Although the model provides an approximate 
reproduction of the continuum in 3C 273 for A < 1450 
A and A > 4000 A, it clearly underestimates the flux 
at 2200 A. In PG 0052-^251, Ton 951, and Mrk 509, 
the models clearly underestimate optical continuum. We 
attribute some fraction of this discrepancy to contamina- 
tion from the host galaxy. This conclusion is supported 
by the apparent anticorrelation between the optical "ex- 
cess" and AGN luminosity as the fractional contamina- 
tion should be weaker in the more luminous AGN if we 
assume an approximately constant host contribution. 

Intrinsic reddening due to dust in the host galaxy may 
contribute to some of the discrepancy between data and 
models in Ton 951 and Mrk 509. However, there is 
disagreement about the form of the reddening curve in 
QSOs, an issue we discuss further in §6.1. As an ex- 
ample, we plot models (dotted curves) which have been 
reddened with the "SMC-likc" reddening curve (Prevot 
et al. 1984) Ax = 1.39 E{B - V){\/ iiia)-'^-^ adopted by 
Richards et al. (2003). We adjust E(B-V) and M in an 
attempt to obtain good agreement at 1450 A and 2200 
A, but still require the flux to match the observed value 
at 4000 A. We find best agreement for E{B - V) - 0.055 
and 0.04 for Ton 951 and Mrk 509, respectively. Al- 
though this improves the agreement in the longward of 
~ 1000 A, the models underestimate the flux at short 
wavelengths. We note that reddening with this extinc- 
tion curve can only worsen the "fit" to 3C 273, since the 
flux at 2200 A is already above the model prediction. 

Table 1 provides a more quantitative comparison using 
spectral slopes. We calculate spectral slopes a {F^, (xu°') 
defined by 

a- (2 \ ^Qg-^^(^'"'^>^) ~^"g-^^(^™") ^ (2) 
V log A 

max 

- log A 

min / 

For the three combinations of continuum windows (1450- 
2200 A, 2200-4000 A, and 1450-4000 A) in Table 1, we 
tabulate three different slopes: one from the data and two 
for models with a* = and 0.9. The a* = slopes cor- 
respond to the models plotted in Figure 1. The a, = 0.9 
models arc selected using the same comparison criterion 
as described above for the a* = models. The a* = 0.9 
slopes are bluer (greater or less negative a) than the 
a* = slopes (cf. Fig. 6), and generally provide a poorer 
match to the observed SEDs. For the best case, PG 
0953-1-414, the the model slopes differ from the contin- 
uum by as much Aa 0.15 for a* = and Aa ~ 0.2 for 
a* = 0.9. For the worst case, Mrk 509, the disagreement 
can be as high as Aa > 1. 

These comparisons suggest that measurement of two 
slopes (between 1450-2200 A and 2200-4000 A) can pro- 
vide a sensible parameterization of the UV continuum. 
They also indicate how problems might arise. Due to 
the small baselines used (A log A = 0.18 and 0.26), a rel- 
atively modest difference in F\ can produce substantial 
discrepancy in the slope. For example, a relative increase 
of 10% in >A at 2200 A would yield A log a = -0.23 
and -F0.16, for the 1450-2200 A and 2200-4000 A, re- 
spectively. This makes the observed slopes particularly 
sensitive to any additional emission which may contam- 
inate the continuum windows. A related concern is that 



a modest amount of dust reddening (cf. Richards et al. 
2003) in the source or host galaxy can substantially alter 
the observed a. For example, E[B — V) =0.03 with 
the above reddening curved leads to Aa = 0.37 and 0.2 
for 1450-2200 A and 2200-4000 A ranges, respectively. 
Therefore, we consider the effects of reddening in §6.1. 

4. SPECTRAL SLOPES OF SDSS QUASARS 
4.1. Sample and Data Reduction 

A principle aim of this work is to look for correlations 
in the continuum of SDSS quasar spectra with M and M. 
For our purposes M can be estimated from the observed 
luminosity, and A'/vir can be inferred from a combina- 
tion of fMg and L3000 discussed in §3.1. These quan- 
tities have already been computed by McLure & Dun- 
lop (2004) for a large sample of sources from the SDSS 
Quasar Catalogue II (Schneider et al. 2003). Therefore, 
we restrict our attention to this sample and refer the 
reader to McLure & Dunlop (2004) for further details. 

From these spectra we must further restrict ourselves 
to two subsamples discriminated by their redshift. At 
longer wavelengths, our slope estimates are calculated us- 
ing the flux measured at A = 2200 and 4000 A. The SDSS 
spectra cover an observed wavelength range from 3800- 
9200 A, so the requirement to observe both wavelengths 
simultaneously restricts us to a range of 0.76 < 2 < 1.26. 
For the shorter wavelength slope, we need coverage from 
A = 1450 — 3000 A in order to calculate the slope and the 
value of isooo which is required for the mass estimate. 
This restricts our second subsample to 1.67 ^ z < 2.09. 
After applying these cuts we are left with 3783 and 2757 
spectra for the low and high z samples, respectively. 

The measurements of fMg and -L3000 in McLure & Dun- 
lop (2004) were performed using spectra from the SDSS 
First Data Release (Abazajian et al. 2003) while our 
spectra were obtained from SDSS Fourth Data Release 
(DR4, Adelman-McCarthy et al. 2006). A recalculation 
of 7;Mg is beyond the scope of this work, biit we have 
recomputed ^3000 with DR4 data, and confirmed that 
the; modest differences in ^3000 have little impact on the 
inferred mass distribution. 

4.2. Spectral Slopes 

In order to compute a, we must extract the flux at 
the wavelengths of interest. We first correct for Galactic 
reddening using the Schlegel et al. (1998) map, and then 
transform to the quasar rest frame. Fluxes and their 
uncertainties are then computed by averaging over 20 
A windows centered on 1450 A, 2200 A, and 4000 A. 
Finally, we use Equation (2) to calculate a from 1450- 
2200 A at high redshift and from 2200-4000 A at low 
redshifts. We use the uncertainties on flux to compute 
the uncertainty in the slope a a, excluding sources with 
Ga > 0.25. This predominantly eliminates sources near 
the flux limit with poor statistics, reducing the samples 
to 3646 and 2706 for low and high redshift, respectively. 
Since these QSOs account for only a small fraction of 
the sample, their exclusion has little effect on the overall 
distribution of a. 

The distributions of a for the two samples are plot- 
ted in Figure 2. The mean slopes for the two samples 
are -0.59 from 1450-2200 A (thin curve) and -0.37 from 
2200-4000 A (thick curve). The differences in the distri- 
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Fig. 2. — Histograms of a calculated for two samples of SDSS 
QSOs. The thin solid curve corresponds the distribution of a mea- 
sured from 1450-2200 A for a sample of QSOs with 1.67 < ^ < 2.09. 
The thick solid curve corresponds the distribution of a measured 
from 2200-4000 A for a sample of QSOs with 0.76 <z< 1.26. 



bution of a between the two samples suggests that the 
QSO spectra may be sUghtly concave with redder slopes 
at shorter wavelengths. This observation is consistent 
with the slope measurement for the relatively nearby 
AGN summarized in Table 1, which are always redder 
at shorter wavc;k;ngths. 

In Figure 3 we plot the mean slope a, binned in terms 
of Mvir and L2200 = ALa evaluated at 2200 A. The error 
bars represent uncertainty estimates for a, correspond- 
ing to a/^/N, where N is the number of spectra which 
contribute to the bin and = J2{'^ ~ ^)^/(-^ ~ 1) is 
the variance in each bin. Differences in the distributions 
of a between the two samples are clearly evident. The 
higher rcdshift. shorter wavelength sample shows weak 
anticorrelations of a with ^2200 ^ind M^ii-. In contrast, 
the lower redshift, longer wavelength sample suggests 
that a is positively correlated with both L2200 and Mvir- 
The strongest trend is the nearly monotonic rise in a as 
^2200 increases. There also appears to be a weaker trend 
in which a increases as Mvir increases to ~ lO^M©, but 
then turns over and begins to decreases at higher masses. 

In order to examine possible correlations between the 
parameters, we plot a 2D distribution of a in the top 
panel of Figure 4. We bin the data in both L2200 and 
Mvir simultaneously. In this plot and several to follow, 
the left and right hand columns correspond to quanti- 
ties measured using the 1450-2200 A and 2200-4000 A 
samples, respectively. In the middle panel we plot the 
standard deviation a, and in the bottom panel we plot N. 
The diagonal dashed lines provide a simple, although rel- 
atively crude, estimate of the Eddington ratio. They are 
curves of constant I/2200/Mvir and would correspond to 
lines of constant Eddington ratio if a linear relationship 
between bolometric luminosity Lboi and ^2200 held for all 
sources. From right-to-lcft each curve represents a factor 
of 10 increase in this ratio and the left-most line would 
correspond to the Eddington limit. For this estimate, we 
assume a ratio iboi/i'2200 = 5, but do not regard this ex- 
act value as particularly significant. Since the mean ratio 
of i22oo/-^3000 ~ 6/5 in our sample, this choice provides 
approximate consistency with the iboi/^aooo = 5.9 rela- 



tion estimated by McLure & Dunlop (2004). 

We find that the 1450-2200 A slopes are not strongly 
correlated with either L2200 or Myir- The highest concen- 
tration of blue slopes are found at low ^2200 and low Mvir 
while the highest concentration of red slopes occurs at 
low to moderate ^2200 and high Mvir (i.e. low L/L^dd)- 
There is a clear trend in the 2200-4000 A slopes in which 
a decreases as L2200 decreases at fixed Myir- A weaker 
variation with Mvir can also be inferred, even at fixed 
luminosity, although it is strongest at high and low Myjr 
where N tends to be lower. At fixed £2200: ct first in- 
creases with increasing Mvir and then decrease at high 
Mvir and low X/I/Edd, consistent with Figure 3. Com- 
parison of a and a in the top two panels shows that the 
typical variance in each bin can be quite large relative to 
the observed trends. The strongest trend in the data, the 
variation in a from L2200 ^ 5 x lO**^ — 10^^ erg s"^ at long 
wavelengths, corresponds to Aa ^ 0.45 while a > 0.3 is 
common. We attribute some, but not all, of this scatter 
to errors in the flux measurements, for which aa ^ 0.15 
is typical. 

Obviously, the variation of a with monochromatic lu- 
minosity depends to some extent on the wavelength used 
to evaluate it. Since it is the ratio of L2200 to ^4000 

determines a in the first place, a positive (negative) cor- 
relation between a and L2200 (^4000) would occur if the 
logarithmic ratio of these luminosities were randomly dis- 
tributed about some mean ratio. We plot the distribution 
of a, binned by ^4000 in Figure 5. As might be expected, 
wc find a weaker, but still positive, correlation between a 
and L4000 than we found between a and L2200 (cf- Figure 
3). 

5. MODEL SLOPES AND MONTE CARLO COMPARISONS 

5.1. Model Slopes 

We now wish to more carefully examine the model 
SEDs described in §2. Although the models we use are 
rather sophisticated in how they treat radiative transfer 
and disk structure, it is useful to begin by considering 
a simple case. Perhaps the simplest spectral model one 
can construct is the multitemperature blackbody. In this 
case one simply calculates the radiative flux in the disk 
as a function of radius and computes a spectrum by in- 
tegrating the emission over the disk surface, assuming a 
Planck spectrum B,^ at the local effective temperature 
Teff. This yields 



U = 477^ j B,{T^s{R))RdR. 



(3) 



To proceed further, one must obtain Teg as a function 
of R. A common approximation is to assume a power 
law form for the flux F oc R~^, yielding 



T'efT = Til 



R 

R\n 



-13/4 



Inserting this form into equation (3) yields 



2AmlaTlh ( hv 



^B^in 



/ 



Xout 2;8//3-l 



- 1 



(4) 



where x = hv/lk-oT). Note that in addition to the explic- 
itly power law dependence on v, Ly is also a function of 
V through the limits of the integral a;i„ = hvj (/cBTln) and 
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Fig. 3. — ID distribution of mean slopes a binned by monochromatic luminosity at 2200 A (left) and viral mass (right). The top 
and bottom panels correspond to slopes measured from 1450-2200 A and 2200-4000 A, respectively. Error bars are calculated assuming 
uncertainties tr/x/iV where c and N axe the standard deviation and number of spectra in eajch bin. 



Xont = hi^ / (kBTout) , For hi^ < ksTout and > fceTln 
we have (x and oc j^'^ exp (—/izy/fceTln), re- 
spectively. At intermediate frequencies fceTln ^ hv ^ 
te^out, the integral is almost independent of v and we 
find oc z^^-s//?, or a = 3 - 8//3. 

For a Newtonian thin disk the flux is given by (Shakura 
& Sunyaev 1973) 

whore R is the radius and Z is a correction factor which 
depends on assumptions about the torque at the inner 
edge of the disk. Typically, T is only a weak function of 
R which approaches unity at large R, so we will ignore 
it for this simple example. Then /? !^ 3 and a ^r^: 1/3 well 
below the peak in the SED. 

For a black hole of mass M, it is useful to scale R 
with Rg = GM/c^ and M with MEdd = -^Edd/c^ = 
47rGM/(cKos), where Kcs is the electron scattering opac- 
ity. With the scahngs m = M/MEdd and r = R/Rg, we 
find 

1/4 

Tin = ( . j:;" , I . (7) 



These results suggest that a ~ 1/3 for hiy <C ksTin, 
with Tin given by equation (7). This result, first derived 
by Lynden-Bell (1969), is commonly referenced as a char- 
acteristic spectral slope for accretion disks (Frank et al. 
1992), but we shall see below that more sophisticated 



models gcnerically give lower values of a for the masses 
and luminosities considered in this work. In part, this 
will be due to the breakdown of the underlying assump- 
tions: that the emission is blackbody and that the flux 
is given by a simple power law with /3 = 3. However, 
it will also not always be the case that hv <^ k^Ti^ at 
the UV wavelengths we are considering. As M increases 
or TO decreases, the frequency at which the SED peaks 
{p ~ SksTin/h) decreases. For hv < fceTin, the spectrum 
should begin to flatten and a is expected to decrease. 

This trend is clearly seen in Figure 6, where we have 
calculated a from 1450-2200 A and 2200-4000 A, using 
equation (2). In the top panel, we compute a for rela- 
tivistic, multitemperature models with a* = 0. We plot 
a as a function L/L^dd (oc m), where L is the model 
bolometric luminosity. The shape of the spectrum, and 
therefore, a depends (weakly) on inclination. Here we 
plot the mean a, averaged over a uniform distribution 
in cosi from cosi = 0.5 to 1. The diagonal dashed and 
dotted curves correspond to lines of constant Tin and 
constant L, respectively. L increases from bottom left 
to top right and Tin increases from the top left to bot- 
tom right. In this case, a can be entirely parameterized 
by Tin- Since the 1450-2200 A slopes are measured with 
continuum windows nearer to the peak in the SED, it is 
always the case that a is lower than measurements made 
from 2200-4000 A. We note that even for these blackbody 
models, a is less than the "canonical" value of 1/3 over 
the entire parameter space of interest. 
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Fig. 4. — 2D distributions binned by monochromatic luminosity L2200 and the viral mass estimate Mvir- The left and right panels 
correspond to the samples with slopes measured from 1450-2200 A and 2200-4000 A, respectively. The top panels display the mean slopes 
o, the middle panels show the standard deviation of a in each bin, and the bottom panels represent the number of spectra contributing to 
each bin. The diagonal dotted lines are approximate estimates for curves of constant Eddington ratio, assuming L^a\ = 5-L2200- See the 
text in 54.2 for further disscussion. 



In the middle panels of Figure 6 we again plot a for 
a, = models, but with spectra based on the atmo- 
sphere models described in §2. For 2200-4000 A slopes, 
we again find that the variations in a can largely param- 
eterized by Tin, although not completely. The primary 
difference is that the slopes are lower than the multitem- 
perature blackbody models and the overall variation with 
Tin tends to be weaker. There are several effects which 
may contribute to these differences, but at these wave- 
lengths the most important modification to the spectrum 
is the Balmer edge at 3650 A. For the majority of annuli 
which contribute to the emission near 4000 A, there IS a 
strong Balmer edge in absorption. The flux emitted by 
any annulus is fixed by the local dissipation of energy, so 
the decrement in flux shortward of the edge is compen- 
sated by an increase longward of the edge. At 2200 A the 
effects of the Balmer edge are minimal. The net effect 
is to increase the ratio 4000 A to 2200 A flux relative to 
what it would be if the edge was absent, yielding a lower 



a. 

The differences in the 1450-2200 A slopes between the 
two models are greater. The contours of constant a are 
generally more vertical. At high M and L/L^dd, they 
actually turn over and nearly follow the lines of constant 
L (dotted curves). Compared with the multitemperature 
blackbody models, we find that a is generally higher, 
particularly at low L, and only at high T is a lower than 
the blackbody prediction. 

The differences between the model slopes at shorter 
wavelengths are largely due to the increasing importance 
of electron scattering opacity in annuli where these pho- 
tons are primarily emitted. As the ratio of scattering 
to absorption opacity increases the spectrum of a partic- 
ular annulus becomes a "modified blackbody" (see e.g. 
Rybicki & Lightman 1979) which peaks at higher pho- 
ton energies. Since total flux is conserved, the increase 
in higher energy photons must be offset by a decrease in 
the flux at lower energies. The emission at any particular 
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wavelength A (in this case 1450 A) is going to eomc from 
a range of annuli with Teg- ~ c/i/(3fcBA). When com- 
pared with blackbody emission, some of these annuli will 
contribute more flux to A and some will contribute less. 
The net eflFect when emission is integrated over all annuli 
will depend on a number of factors: how close Tcff is to 
Tin, the wavelength at which k^s becomes the dominant 
opacity, and how strongly the spectrum is deformed from 
blackbody. For low values of Tin and low L this deforma- 
tion tends to enhance the flux at 1450 A, but has little 
effect on the flux at 2200 A, increasing a. 

In the bottom panels, we again plot slopes for atmo- 
sphere based spectra, but with a* = 0.9 instead of a* = 0. 
The dependence of a on L/L-E^d and M is very similar to 
the o* = case in that contours of constant a have sim- 
ilar shapes. However, the spectra arc gcncrically bluer 
since a is larger than the a* = slopes over the whole 
range covered by the plot. There are multiple ways in 
which changing a* modifies a, but the dominant effect 
is the shift in the peak of the SED. In these models, the 
inner radius r^n is determined by the location of the in- 
nermost stable circular orbit (ISCO), which is smaller for 
larger a*. From equation (7), we can see that this makes 
Tin larger. At fixed L/L-Edd, "rn must also decrease to 
offset the increase in radiative efficiency, which would re- 
duce Tin- However, the former effect dominates, and the 
peak of the SED shifts to shorter wavelengths, generally 
increasing a. 

5.2. Monte Carlo Comparisons 

With accurate mass, spin, and inclination estimates, 
we could compare the slopes in Figure 4 directly with 
the models, such as those shown in Figure 6. However, 

we do not have any practical means of estimating either 
the spin or the inclination. Furthermore, we do not know 
with ('Certainty the; accuracy of mass estimates M^{^ ob- 
tained with equation (1). In order to account for these 
uncertainties we construct slope distributions similar to 
those shown in Figure 4 using Monte Carlo methods. 

In order to characterize the distributions of model pa- 
rameters we must make several assumptions. First we 



must adopt a distribution of inclinations. In absence of 
obscuration, a uniform distribution in cos i would be ex- 
pected from isotropic emission. However, the disk models 
do not produce isotropic emission. Disks viewed nearly 
face on (cosi ~ 1) have larger fiuxes than edge-on disks, 
and should be somewhat enhanced near the flux limit in 
a flux limited sample. Furthermore, angle dependent ob- 
scuration is a fundamental tenet of the unifled model of 
AGN (Antonucci 1993). In the unified model, obscuring 
material lies in the disk plane (the "torus"). Broad emis- 
sion line objects, such as the type I QSOs discussed here, 
will be viewed nearly face on, up to the opening angle of 
the torus. Here, we adopt a uniform distribution in cos i 
from cosi = 0.5 to 1, wlicre tlie lower limit corresponds 
to an opening angle of 60°. Since the differences in cosi 
produce at most a factor of two in flux, the assumption 
of a uniform distribution will only produce a small error 
near the flux limit. 

We must also specify a*. Unfortunately the distri- 
bution of a* is highly uncertain. Estimates for a* in 
AGN are limited to a handful of bright, relatively nearby 
sources with broad Fe Ka lines. In some cases, such as 
MCG -6-30-15, very high values of a* (~ 1) are inferred 
(see e.g. Tanaka et al. 1995). A combination of empir- 
ical and theoretical arguments seem to favor a* --^ 0.9 
(Gammie et al. 2004), although there are considerable 
uncertainties. Given the large uncertainties, we simply 
choose two characteristic spins, a, = and 0.9, as repre- 
sentative examples. 

Finally, we must account for the uncertainties in Mvir- 
Although these uncertainties are diflicult to estimate ro- 
bustly, some insight may be obtained by comparing es- 
timates which utilize the H/3 line width and monochro- 
matic luminosity at 5100 A. McLure & Dunlop (2004) 
compute the ratio log(Mvir[H/3]/Mvir[Mg II]), flnding a 
dispersion of 0.33 dex. This, of course, does not ac- 
count for any systematic errors common to the two meth- 
ods. The Mvir [H/3] estimates rely on luminosity based 
-Rblr estimates which are calibrated by matching rever- 
beration mapping estimates. For example, Vestergaard 
(2002) find that 70% of Mvir [H/3] mass estimates match 
the reverberation mapping estimates to within a factor 
of 3. The reverberation mapping estimates, in turn, are 
claimed to have a typical precision of ^ 30% (Peterson 
et al. 2004). Here, we adopt 0.4 dex as a fiducial value of 
the uncertainty in log Mvir. Given the scatter in the re- 
lations discussed above, we view this as a lower limit on 
the typical error. Therefore, we also consider the impact 
of assuming a larger miccrtainty (0.8 dex) below. 

With these assumptions we can begin computing 
Monte Carlo slope distributions. We start by assum- 
ing a distributions of T2200, ^Mg, and Mvir identical to 
those in §4.2. Next, we use random deviates to draw val- 
ues of the "actual mass" Mo and cosi. Here, Mo is used 
to account for possible errors in the mass estimates. It 
is drawn from a log normal distribution with a mean of 
Mvir and gm = 0.4 dex. This prescription means that 
the mass distribution of the models will be somewhat 
broader than the distribution of A'/vir- However, there is 
some evidence for a real cutoff in the mass distribution 
at high mass. Therefore, we enforce a maximum mass 
logMo/M© < 9.5. Then, for each choice of Mq and i, a 
value of L/TEdd is chosen so the model monochromatic 
luminosity matches i2200- Finally, given these values 
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constant Tin a^nd the diagonal dotted lines represent lines of constant L. See the text in §5.1 for further discussion. 



of L/LEdd, Mq, i, and a choice of a* we can compute 
a for the corresponding model. These slopes are then 
binned in exactly the same manner as the data, allowing 
us to compare with the observed a and a. (Note that 
the distributions of TV are identical to the observations 
by construction.) 

We plot the resulting 2-D distribution of a for — 
in the top panels of Figure 7. For both long and short 
wavelength slopes, there are clear mass and luminosity 
dependent trends in a. At fixed L220O1 ct decreases with 
increasing Afvir and at fixed Mvir, a generally decreases 
with decreasing ^2200- Since L2200 is roughly propor- 
tional to the bolometric luminosity, the trends are as 
expected from middle left panel of Figure 6. Compari- 
son with Figure 4 shows that these strong variations are 
clearly discrepant with the weak trends in a inferred from 
the observations. Although we do not plot the results, 
we have performed the equivalent exercise using the mul- 



titemperature blackbody slopes. The model trends are 
even stronger, due the rapid reddening of the spectra at 
high masses and low luminosities. The discrepancies are 
particularly large at short wavelengths, where the mul- 
titemperature blackbody slopes are significantly redder 
than the observed spectra. 

In the the top panel of Figure 8, we plot the distribu- 
tion of standard deviations corresponding to the slopes in 
Figure 7. The distributions of a are largely determined 
by the assumed uncertainty in Mvir ■ If the assumed 
uncertainty ctm was identically zero, each a would be al- 
most completely determined by L2200 and Mvir since the 
variation in a with cosi is relatively weak. In that case, 
cr < 0.05 would be typical. However, comparison with 
the middle panel of Figure 4 shows that despite this scat- 
ter induced by the mass uncertainties, a for the model 
slopes is still substantially below the standard deviations 
obtained from the observed slopes. 
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Fig. 7. — Simulated 2D distributions of a binned by L2200 and Mvir for comparison with the top panels of Figure 4. The left and right 
panels correspond to the slopes measured from 1450-2200 A and 2200-4000 A, respectively. The distributions of a are calculated from thin 
disk models using a Monte Carlo method described in §5.2 to approximate the distribution of mass and luminosity. Each row corresponds 
to a different set of assumptions for the black hole spin and the uncertainties in the mass estimates. In the top and middle panels we plot 
a calculated for models with a« = while the bottom panels show the results for a, = 0.9. The top and bottom panels assume aj^i = 0.4 
dex, while the middle panels assume (tm = 0.8 dex. 



Given the nature of the discrepancy between models 
and data, an obvious concern is the effect of errors in 
the mass estimates. If the errors in the mass estimates 
are larger, they might smear out any mass dependence in 
the observed slopes and increase a. Of course, the actual 
form of the mass error distribution remains an important 
source of uncertainty, and the log normal distribution 
employed here may not be an adequate approximation. 
However, since we have no reliable means of indepen- 
dently measuring this distribution, we simply accept this 
as a caveat and consider the effects of larger errors by in- 
creasing (Tj\/ from 0.4 to 0.8 dex. We plot the results in 
the middle panels of Figures 7 and 8. 

In the middle panel of Figure 7, we find that the re- 
sulting distributions of a are still incapable of repro- 
ducing the observations, but increasing gm to 0.8 dex 
does smear out the strong variation seen the ctm = 0.4 
dex expected. At short wavelengths, this re- 

moves almost all of the variation in the slopes with Mvir 



and -^2000 1 better agreement with the lack of a strong 
trend in the observed slopes. However, the mean slopes 
are still larger than those that are observed. At longer 
wavelengths, there is still some trend for larger slopes for 
higher L/ L-em, which is inconsistent with the observed 
distribution. For L/L-^dd ^ 0.1, the model slopes are too 
large at low L2000 smd too low at high £2000- 

In the middle panels of Figure 8, we find that a also 
increases, as expected. This provides better agreement 
with observed a, although a is still generally too low. We 
also find that a decreases as L/L-^dd increases, a trend 
not seen in the data. This decrease in a is also consistent 
with Figure 6, which shows that a is a weaker function of 
L/L-Edd and M for the low M and high L/L-^dd models 
which predominantly contribute to the low a bins. 

Except at high luminosities and long wavelengths, we 
nearly always find that the models predict a larger 
(bluer) than the observed values. A comparison of the 
bottom and middle panels of Figure 6 suggests that in- 
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Fig. 8. — The standard deviation a for the simulated distributions of a plotted in Figure 7. Each row correspond to a different set of 
assumptions for the black hole spin and the uncertainties in the mass estimates, as described in the caption of Figure 7. 



creasing a» will only make the discrepancy worse, since 
the a* = 0.9 models slopes are everywhere larger than 
the o, = case at equivalent M and L/L-Edd- Indeed, 
this is precisely what we find in the bottom panels of 
Figures 7 when we plot a for a, — 0.9 and ctm = 0.4 
dex. For the range of Mvir and ^2000 sampled, a varies 
less strongly than in the — Q case. However, this 
also means that a is even lower, increasing the disagree- 
ment with the data. Thus, it appears that the a, = 0.9 
models are an even poorer match to the data that the 
a* = models. This result may be problematic since, as 
discussed above, there is some evidence which suggests 
a, > 0.9 may be more common than a, ~ 0. 

6. DISCUSSION AND CONCLUSIONS 

6.1. The Effects of Dust Reddening on Spectral Slopes 

It is clear that some additional process (or processes) 
must be modifying the spectrum if the models are to be 
reconciled with the observed slopes and we discus a num- 
ber of possibilities below. One of the main discrepancies 
is that the observed spectra are generally redder than 
the models predict. One likely possibility is that the 



SEDs are altered by wavelength dependent extinction 
from dust intrinsic to the source or host galaxy. Dust 
emission can be seen in the infrared in many AGNs, and 
is almost certainly present at some level in our sample 
(see e.g. Richards et al. 2003). The uncertainties are 
mainly questions of how much dust is present and what 
is the wavelength dependence of the extinction. 

We can incorporate the effects of dust reddening di- 
rectly into our Monte Carlo distributions by simulating 
the effect of extinction on the model SEDs. Multiplying 
our artificial spectra by a reddening curve will change 
both the values of a and £2200 inferred by an observer. 
In order to proceed we must specify the shape of the 
reddening curve and specify the amount of extinction. 
Often, the degree of reddening is parameterized by the 
color excess E(B-V) which is the difference in extinction 
between the B and V bands, expressed in magnitudes. 
Since this distribution is unknown, we simply adopt a 
uniform distribution of color excess between zero and 
some specified maximum value. 

The wavelength dependence of the extinction is also 
uncertain. The determination of QSO reddening curves 
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is an area of active research and ongoing debate 
(Richards et al. 2003; Gaskell et al. 2004; Hopkins et al. 
2004; Czerny et al. 2004; Willott 2005). Much of the 
difficulty in determining these reddening curves stems 
from the problem of disentangling intrinsic variations in 
the quasar SED (lines or continuum) from the effects of 
dust. This makes it difficult to unambiguously identify 
"unreddened" AGN sources for comparison with a "red- 
dened" population. Despite this uncertainty, it is widely 
accepted that the ~ 2200 A feature seen in the Galac- 
tic interstellar medium (see e.g. Cardelli et al. 1989) is 
extremely weak or absent in AGN reddening curves. 

Disagreement primarily arises over the steepness of the 
reddening curve, particularly at wavelengths shortward 
of ~ 3000 A. For example, GaskeU et al. (2004) derive a 
flat far UV reddening curve, while the analysis of Hop- 
kins et al. (2004) favors a steeper, SMC-like curve. An- 
other calculation by Czerny et al. (2004), which is based 
on ratios of the color selected SDSS quasar composites 



of Richards et al. (2003) , prefers a curve which is in be- 
tween: somewhat flatter than the SMC, but still steeper 
than that found by Gaskell et al. (2004) . 

Clearly, if extinction is important, the observed slope 
distributions will depend significantly on the form of the 
reddening curve. Therefore, we have considered three 
possibilities: the approximate forms of the curves de- 
rived by Gaskell et al. (2004) (see their Appendix) and 
Czerny et al. (2004) (equation [3]), and the SMC-like 
curve used by Richards et al. (2003) (^a = 1-39 E[B - 
y)(A//im)-i-2). 

In Figure 9 we plot Monte Carlo distributions equiv- 
alent to those in Figure 7, but including the effects of 
dust. For this example we use the SMC-like reddening 
curve and choose maximum E(B-V) so that short wave- 
length slopes approximately match the observed slopes 
(cf. Figure 4). This matching requires E(B-V) ^ 0.03 
and 0.055 for a* = and 0.9, respectively. These mod- 
est values would not contradict the results of Richards 
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et al. (2003) who found that 94% of their sample would 
be consistent with E(B-V) < 0.04 for the same redden- 
ing curve. We note that the modification of L2200 due 
to extinction is typically small for the required range of 
E(B-V). In principle, a larger amount of extinction could 
introduce a correlation where dust reddened sources are 
less luminous than bluer, unreddened spectra, and pos- 
sibly explain the observed trend of higher a for higher 
^2200- Also, the small amount of extinction at 2200 A 
introduces considerably less scatter in a than the uncer- 
tainties in the mass estimate. Therefore, we do not replot 
the standard deviation of the bins, since the plots differ 
little from those of Figure 8. 

The reddening curve of Czerny et al. (2004) has a flat- 
ter wavelength dependence than the SMC-like curve, and 
therefore requires greater extinction to produce the same 
amount of reddening. In order to match the short wave- 
lengths slopes of Figure 4, we require a maximum E(B-V) 
~ 0.6 and 0.12 for a* = and 0.9, respectively. The rela- 
tive flatness of this reddening curve in comparison to the 
SMC-like curve also leads to a greater degree reddening 
at longer wavelengths to produce the same amount of 
reddening a short wavelengths. This is particularly sig- 
nificant for a, — 0.9 model. Since the observed slopes 
are redder at lower ^2200, this improves the agreement 
with the data for low L2200 ) but provides a poorer match 
for higher values of L2200J where slopes arc bluer. 

The Gaskell et al. (2004) curve is nearly flat shortward 
of 3000 A. Therefore, it requires substantially greater 
extinction to compensate for the slope discrepancy be- 
tween the models and data at short wavelengths. We 
cannot easily parameterize the amount of extinction re- 
quired because it is large enough to boost the intrinsic 
(unreddened) luminosity by an order of magnitude or 
more. This would make the models substantially super- 
Eddington for the assumed masses, exceeding the upper 
limit of our grid where the thin disk assumptions are no 
longer self-consistent. 

Although our assumed distribution of color excess 
which is uncorrelated with luminosity can improve the 
agreement between the models and data at short wave- 
lengths, it cannot reproduce the L2200 dependent trend 
at longer wavelengths. To some degree, this conclusion 
hinges on our assumption that the amount of dust ex- 
tinction is independent of luminosity. If, for example, 
the amount of dust extinction anticorrelates with lumi- 
nosity, the most luminous QSOs could have unreddened, 
intrinsically blue spectra while lower luminosity QSOs 
would have lower slopes due to the dust reddening, in 
agreement with the observed trend. However, the lack 
of a similar luminosity dependent correlation at short 
wavelengths would then require a reddening curve which 
is flat shortward of ~ 2200 A, such as that proposed by 
Gaskell et al. (2004). Since the short wavelength slopes 
are redder than predicted by the models, matching the 
data with the models at short wavelengths is not possi- 
ble with such a flat reddening curve. Therefore, even if 
such luminosity dependent redding is plausible, it seems 
unable to simultaneously account for both the short and 
long wavelength spectral slopes if the underlying con- 
tinuum is well approximated by the models. However, 
such a scenario may account for the observed trend if 
the short wavelength continuum is intrinsically redder 



than the models predict. 

6.2. The Effects of Irradiation 

Self-irradiation presents another possibility for explain- 
ing the red slopes in these systems. In fact, correlated 
variability on timescales comparable to the light travel 
times suggests that irradiation must be occurring at some 
level (see e.g. Krolik et al. 1991). In the UV, we expect ir- 
radiation to become increasingly less important for deter- 
mining the spectrum as we move to shorter wavelengths 
in the UV, since local dissipation must dominate the flux 
near the peak of the SED. Nevertheless, irradiation may 
still play some role so wc briefly consider its effects by 
examining a simple model. 

The slope modification due to reddening will be 
strongly dependent on the geometry of the accretion fiow. 
A simple estimate of the slopes of irradiated disks may be 
obtained from equation (5) which implies a = 3 — 8//3 for 
hv <C fceTin. In order to estimate /?, we need to specify 
the radial dependence of the irradiated flux Fi^ oc R^^ . 
With simple geometric arguments, one can show (see e.g. 
Blaes 2004b) that 



_ L41 - a) fH\ f dlogH _.,HA 
AnB? \r) \d\ogR RJ' 



(8) 



In order to obtain this relation, wc approximate the irra- 
diating continuum (the inner disk, or possibly a corona) 
as a point source at i? = and a height above the 
midplane. The disk surface at the point of irradiation 
is parameterized by the height H{R) at radius R. The 
albedo a may also be a function of R, but we will ignore 
this dependence for simplicity. 

For iJ» > H{R), we find Fi„ a R'^. This depen- 
dence is easily understood by noting that emission from 
the point source falls off as and is intercepts the 
disk surface with angle such that cos 6* = H^/R. Since 
H^, is independent of R, we find F^^^ oc R~^ as inferred 
from equation (8). If this reprocessed fiux dominates the 
locally dissipated flux and is reradiated as a blackbody 
we again have /3 = 3 and a = 1/3, equivalent to the bare 
thin disk case. 

If ff, < H{R), the disk must flare in order to produce 
significant irradiation. If we parameterize H as a power 
law {H (X R^), we find /3 = 3-ror a = (1 -3r)/(3-r). 
An important case is F = 1 for which /3 = 2 and a = —1. 
For F > 1, a < — 1 and vice- versa. This suggests that in 
the region of parameter space where most of our slopes 
lie, — 1 < a < 1/3, we require F < 1. Such models are 
concave in shape and self-shielding at sufficiently large 
radii. 

Obtaining slopes in this range via irradiation may 
therefore require some level of fine tuning. One possi- 
bility is that the disk transitions from a flat or convex 
(F > 1) solution to a concave (F < 1) solution at the 
range of radii which give rise to the UV emission so that 
self-shielding occurs only at larger radii and longer wave- 
lengths. A second possibility is that there is flaring with 
F > 1, but that the reprocessed flux does not dominate 
the local dissipation. The differing fraction and radial 
dependences of the local dissipation and reprocessed flux 
give rise to a range of a between ~ — 1 and ~ 1/3. Such 
a concurrence would be somewhat surprising, because a 
comparable contribution from the locally dissipated and 
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reprocessed emission will only occur over a limited range 
of radii, due to their different radial dependences. 

Nevertheless, self-shielding geometries may be useful 
for explaining some aspects of the observed slope dis- 
tribution. In the top panel of Figure 4 we see that for 
the highest luminosity sources, slopes are, on average, 
bluer at longer wavelengths than at short wavelengths. 
The redder than expected slopes at short wavelengths 
could be produced by irradiation of a portion of the disk 
surface which blocks emission from reaching larger radii. 
The unirradiated flow at larger radii and lower Tes could 
remain dominated by the local dissipative flux and pro- 
duce an intrinsically bluer slope. 

In order for such interpretations to be valid, the ge- 
ometry of the irradiated accretion flows needs to be ex- 
plained. Our thin disk models for supermassive black 
holes arc radiation pressure dominated in their inner- 
most radii (Shakura & Sunyaev 1973), implying a scale 
height which is nearly independent of radius. With this 
model, we would not expect reprocessing to significantly 
alter the spectrum until the disk transitions to the gas 
pressure dominated regime where H cx i?2i/20 (shakura 
& Sunyaev 1973), yielding a ^ — 1 if reprocessing domi- 
nates the local flux. The transition radius from radiation 
to gas pressure dominance depends on M, ass, '>^, and 
a*, but is generally located at 200-400 Rg. 

Since the bulk of the UV radiation in our models is 
radiated at radii < 200i?g, some other mechanism must 
be modifying the disk structure in order to explain the 
red slopes via irradiation. One possibility is that the ver- 
tical extent of these disks may be substantially modifled 
due to the magnetic support, as suggested by numeri- 
cal simulations (see e.g. Turner 2004; Hirose et al. 2006). 
However, at present, such calculations remain too un- 
certain to yield a predictive model for the reprocessing. 
Another possibility is that backscattered radiation from 
an outflow might modify the spectnmi (e.g. Nikolajuk 
et al. 2004), but this will depend on the (unknown) out- 
flow geometry, even if such outflows prove to be common. 

6.3. The Dependence of Slope on Luminosity 

As discussed in §5.2 and §6.1, the luminosity depen- 
dent slopes at long wavelengths present a challenge for 
the Hubeny et al. (2000) models since the discrepancies 
cannot be simply attributed dust reddening or errors in 
the mass estimates. As discussed in §6.1, adding an ad- 
hoc luminosity dependence to the dust extinction might 
account for the dependence of a on L2200 at long wave- 
lengths, but not without creating problems at shorter 
wavelengths. 

Due to the short baselines involved (log(4000/2200) = 
0.26 and log(2200/1450) = 0.18) the slopes can be sub- 
stantially modified by relatively small changes in fiux of a 
continuum window. Therefore, we consider "contamina- 
tion" from the host galaxy and/or the BLR emission as a 
possible solution. In such a scenario, the observed trend 
would imply either an increasing contribution to fiux at 
4000 A as ^2200 decreases or an increasing contribution 
to the flux at 2200 A as L2200 increases. Inspection of 
the bottom, left panel of Figure 3 shows a variation in 
the mean slope Aa ~ 0.45 from ^2200 < 10^^ erg s~^ to 
■^^2200 ^ 10^^ erg s~-^. This is only shghtly greater than 
the typical standard deviation in the individual bins and 
corresponds to a 30% variation in the relative flux 



between the two continuum windows. 

One scenario would be a contribution from the host 
galaxy which is roughly independent of the QSO lu- 
minosity and contributes 30% at low Z/2200- Based 
upon the strength of absorption lines in their compos- 
ite spectrum Vanden Berk et al. (2001) estimate a 7%- 
15% contribution at the location of Ca II A3933 and Na I 
A5896. Therefore, the trend cannot be explained entirely 
by host galaxy contamination. A second possibility is 
excess BLR emission at 2200 A which correlates with 
Z/2200- Inspection of the Vanden Berk et al. (2001) com- 
posite spectrum indicates emission from high excitation 
Fe II lines is probably the main contaminant at ~ 2200 
A. However, a comparison of composite spectra for QSOs 
with 1^2200 above and below 3 x 10'''^ ergs~^ does not pro- 
vide evidence for such a large change in the equivalent 
width of these lines so we conclude that this explanation 
is also unlikely. 

It is also possible that the underlying accretion flow 
is thin disk, but that the Hubeny et al. (2000) models 
have not properly accounted for the emission near ~ 4000 
A. The model grid only extends to Tdr = 10^ K, so 
the emission for annuli in the disk at larger radius and 
lower Toff are approximated by blackbodies. For models 
with Tes < 10^ K density inversions can occur due to a 
fall in electron pressure when H recombines (see Fig. 9 
of Hubeny et al. 2000). This makes it very difficult to 
predict the true equilibrium structure, and therefore the 
spectrum, of an annulus in a real, turbulent accretion 
flow. 

The lack of models Tefj < 10* K creates a problem 
because these annuli may still have significant Balmer 
edges. The presence of edge at ^ 3650 A tends to create 
an excess of emission longward of the edge and a deficit at 
shorter wavelengths relative to the predictions of a pure 
blackbody. As can be seen by Figure 11 of Hubeny et al. 
(2000), replacing non-LTE model spectra with blackbod- 
ies may produce a reduction in the fiux at 4000 A. This 
result, coupled with the observation that the strength of 
the edge is anticorrelated with luminosity at fixed mass 
(see Figure 13 of Hubeny et al. 2000), suggests that the 
Balmer edge may play a greater role in producing the 
observed trend than the current models predict. A ro- 
bust determination of the spectrum at these wavelengths 
might ultimately require spectral models coupled with 
realistic disk simulations, and is beyond the scope of this 
work. 

It is also possible that selection biases and incomplete- 
ness play some role in producing the observed trends. 
The spectroscopic targeting of QSOs in SDSS is discussed 
in Richards et al. (2002). For this sample, we consider in- 
completeness near the flux limit to be the greatest cause 
for concern since the majority of our sources are well 
separated from the stellar locus. 

The spectroscopic targeting algorithm rejects sources 
based upon i band magnitude. For ugri and griz selected 
sources the magnitude limits are i < 19.1 and i < 20.2, 
respectively. The i band filter is approximately centered 
at ~ 7600 A, which corresponds to 3800 Aand ^ 2500 
Afor the low (2 ~ 1) and high (z < 2) redshift samples 
respectively. This could lead to a deficit of blue quasars 
at low £2200 in the long wavelength, low redshift sample. 
For an equivalent L220O) red quasars will have a larger 
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flux near 4000 A than sources with bluer slopes, and 
would be more likely to make it into the spectroscopic 
sample if they are near the i flux limit. We have consid- 
ered this possibility by plotting the slopes as function of 
i for our sample, but we do not find a significant deficit 
of blue slopes near the flux limits. In fact, we find a de- 
creasing fraction of red quasars for i < 17.5. This deficit 
of red spectra with large fiux in i seems to be due to a 
real paucity of red sources with high luminosities. 

6.4. Comparison with Previous Work 

There have been a number of previous observational 
tests of the Hubeny et al. (2000) models considered here. 
For the most part, these comparisons have involved in- 
dividual sources or a relatively small sample with broad- 
band spectral coverage. Blacs et al. (2001) fit models to 
spectra of 3C 273, finding poor agreement in the near 
UV. The reverberation mapping mass estimate has since 
been revised, bringing the model and data into better, 
though not perfect agreement (cf. §3.2). 

Blaes (2004a) calculated spectral slopes as a function 
of mass and Eddington ratio for a subset of sources in the 
Shang et al. (2005) sample. The results showed signifi- 
cant scatter and rather poor agreement with the model 
predictions, but the analysis only included a few dozen 
objects. Using the same data, Shang et al. (2005) com- 
pared the optical, optical-to-UV, and far UV spectral 
slopes to the models, finding rough agreement. A recent 
analysis of SDSS quasars observed by the Galaxy Evo- 
lution Explorer (GALEX) (TrammcU et al. 2007) also 
examined the far UV properties of AGNs, again finding 
evidence for slope changes near ~ 1000 A, m approxi- 
mate agreement with the model predictions. 

A recent analysis by Bonning ot al. (2006) relates most 
directly to our current work. They use another, overlap- 
ping sample of SDSS QSO spectra and measure color 
ratios between continuum windows located at 1350 A, 
2200 A, 4000 A, and 5100 A. They also compare their 
results with the a selection of Hubeny et al. (2000) mod- 
els chosen to approximate the distribution of Eddington 
ratios and masses inferred from the data. 

The choice of variables used for binning is one of the 
major differences between the Bonning et al. (2006) anal- 
ysis and this work. Bonning et al. (2006) consider the 
evolution of observed and artificial spectra as a function 
of a single variable corresponding to Tin. For the vi- 
ral mass estimates used here, the dependence of Tin on 
VMg and 1^3000 can be obtained by combining equations 
(1) and (7). If we make the additional assumption that 
Lboi oc L3000, we find 

-lin oc -1^3000 '^Mg- 

Following McLure & Dunlop (2004), we have used S = 
0.62 which yields a weak dependence on luminosity. Bon- 
ning et al. (2006) use S = 0.5 exactly so that Tin is a 
function of VMg alone. 

To facilitate comparison, wc have rcplottcd our dis- 
tributions as a function of Tin in Figure 10, using flux 
ratios in place of a on the vertical axis. Following equa- 
tion (5) of Bonning et al. (2006), we use the relation 
logTin = 5.43 — log(wMg/3000 kms~-^) for the horizon- 
tal axis. We plot the distributions of both the observed 
(thin curve) and model (thick curve) flux ratios. The 



model curves are calculated using the same Monte Carlo 
distribution as in the top panel of Figure 7. At long 
wavelengths the observed SEDs are reddest at low and 
high Tin. This differs from the model which are red at 
low Tin and become monotonically bluer as Tin increases. 
At short wavelengths the observed flux ratios are roughly 
independent of Tin, while the model fluxes are again red- 
dest at low Tin and become monotonically bluer as Tin 
increases. These results are qualitatively consistent with 
Figure 3 of Bonning et al. (2006), who measure colors at 
1350 A as opposed to the 1450 A window used here. 

Bonning et al. (2006) infer that the observed reddening 
at high Tin may be related to Eddington ratio since most 
of the objects contributing to the highest Tin bins have 
L/L-Edd ^ 0.3. This is also roughly consistent with our 
findings. In our sample this result can be understood by 
examining the top loft panel of Figure 4 and considering 
equation (7). Figure 4 shows that most of the QSOs in 
our sample do not radiate significantly above the Edding- 
ton limit. Equation (7) implies that the highest values of 
Tin are obtained for low values of M and high Eddington 
ratios (i.e. high m). As a result, the sources in our sam- 
ple with the highest Tin (and, therefore, low vyig) tend to 
occupy the low Mvir, low L2200 corner of the plot. Since 
the observed slopes are predominant functions of T2200 
with redder slopes at lower luminosities, these QSOs also 
tend to be redder than average. Of course, bins with 
slight lower Tin also include low T2200 objects, but the 
average a is still larger due to the increasing fraction of 
higher Myir and T2200 sources, which tend to be bluer. 
Therefore, the problem of understanding the redding at 
low Tin is intimately connected to the question of why 
the mean slopes are predominantly functions of T2200 
(rather than T220o/-^vir) at long wavelengths which we 
discussed in §6.3: 

6.5. Conclusions 

We have shown how slopes from artificial SEDs of thin 
accretion disks (Hubeny et al. 2000) vary with black hole 
mass and bolometric disk luminosity. As expected from 
naive models, we find that the slope a generally decreases 
as M increases at fixed L. We have shown that the UV 
spectral slopes of models based on radiative transfer cal- 
culations differ measurably from those with simple black- 
bodies, and considered how the slopes are modified by 
changes in the black hole spin. 

We first compared the models against five broadband 
SEDs of nearby, bright AGN (Shang et al. 2005) which 
also have reverberation mapping mass estimates. For 
the more luminous sources, the models can roughly re- 
produce the observed flux ratios in continuum windows 
at 1450 A, 2200 A, and 4000 A. At lower luminosities, the 
short wavelength slopes are substantially redder than the 
model predictions, and may indicate significant redden- 
ing by dust local to the source or host galaxy. We then 
measured a for 6352 QSOs, using these same continuum 
windows. We find only a weak trend with mass when 
virial estimates are used at short (1450-2200 A) or long 
(2200-4000 A) wavelengths. Even if we allow for errors 
in the mass estimates with a log normal distribution and 
a standard deviation of 0.4 dex, a much stronger mass 
dependent trend is observed in the model slope distribu- 
tions which is not consistent with the observations. 
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Fig. 10. — Distribution of flux ratios binned as a function of Tin for comparison with the results of Bonning et al. (2006). The right 
and left panels correspond to the slopes measured from 1450-2200 A and 2200-4000 A, respectively. We calculate Tin from v^g using the 
relation logTi^ = 5.43 — log(t)Mg/3000 km s^-*^). The top axis displays the value of logv^g for comparison with Figure 3. We plot the 
distribution of the observed SDSS flux ratios as a thin curve. The solid curves show distributions of model flux ratios based on the Monte 
Carlo distributions. The plotted curves correspond to the same distribution plotted as a 2D distribution in the top panel of Figure 7 with 
a* = and <tm = 0.4 dex. 



A possible explanation is that the mass estimates typ- 
ically have larger errors than our assumed distribution 
predicts. In support of this possibility, we find that in- 
creasing the typical mass estimate error to ^ 0.8 dex is 
sufficient to erase most of the mass dependent trend in 
the model slopes. However, this increase alone is insuffi- 
cient to bring the slope distribution into agreement with 
the data. 

The multitemperature blackbody models yield slopes 
which are too red at short wavelengths for black holes 
with M > IO^Mq. With the exception of the longer 
wavelength slopes of the most luminous QSOs, we al- 
ways find that the observed slopes are redder than 
the non-LTE atmosphere based model predictions for 
Schwarzschild black holes. The discrepancy is even 
greater for spinning black holes which are generally bluer 
than their non-spinning counterparts. We suggest that 
much of this discrepancy could be accounted for by dust 
reddening in the source or host galaxy. If we use an SMC- 
like reddening curve (Richards et al. 2003), we require 
E(B-V) < 0.03 and 0.055 in order to obtain agreement 
between the short wavelength slopes. Reddening curves 
which arc flatter at wavelengths shortward of 2200 A (e.g. 
Czcrny et al. 2004; Gaskell et al. 2004) would require a 
greater color excess. If the reddening curve is as flat as 
that derived by Gaskell et al. (2004), dust reddening will 
have little impact on the short wavelength slopes, and 
the observed slopes should very nearly match the intrin- 
sic slopes. In that case, the models considered here would 
not be consistent with the observed slopes. 

At longer wavelengths (2200-4000 A), the observed 
slopes are generally bluer at high L2200 and redder at 
lower luminosities, for a fixed ratio relative to the Ed- 
dington luminosity. The models, however, predict a 



trend which is predominantly determined by Edding- 
ton ratio and arc not consistent with this result. This 
discrepancy remains even after we simulate the effects 
of dust reddening. This discrepancy may partly arise 
from difficulties in properly modelling emission near the 
Balmer edge which can affect the flux at 4000 A sig- 
nificantly. Improving the models at these (and longer) 
wavelengths is an important step for future studies. 

Overall, we find no clear signature of bare, thin accre- 
tion disks from the distribution of observed UV spectral 
slopes. Nevertheless, we do not believe the present anal- 
ysis is sufficient to rule out a dominant contribution from 
such models. We consider uncertainties in the amount of 
extinction, wavelength dependence of reddening curve, 
and precision of the mass estimates to be the most im- 
portant caveats. If the mass estimates prove to be suf- 
ficiently precise (i.e. correct to within a factor of four) 
or dust reddening proves sufficiently weak, the actual ac- 
cretion flows must differ significantly from the models 
employed here. 
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